Heat System¶
District heating with thermal storage and time-varying prices.
This notebook introduces:
- Storage: Thermal buffer tanks with charging/discharging
- Time series data: Using real demand profiles
- Multiple components: Combining boiler, storage, and loads
- Result visualization: Heatmaps, balance plots, and charge states
Setup¶
In [1]:
Copied!
import numpy as np
import pandas as pd
import plotly.express as px
import xarray as xr
import flixopt as fx
fx.CONFIG.notebook()
import numpy as np import pandas as pd import plotly.express as px import xarray as xr import flixopt as fx fx.CONFIG.notebook()
Out[1]:
flixopt.config.CONFIG
Define Time Horizon and Demand¶
We model one week with hourly resolution. The office has typical weekday patterns:
In [2]:
Copied!
# One week, hourly resolution
timesteps = pd.date_range('2024-01-15', periods=168, freq='h')
# Create realistic office heat demand pattern
hours = np.arange(168)
hour_of_day = hours % 24
day_of_week = (hours // 24) % 7
# Base demand pattern (kW)
base_demand = np.where(
(hour_of_day >= 7) & (hour_of_day <= 18), # Office hours
80, # Daytime
30, # Night setback
)
# Reduce on weekends (days 5, 6)
weekend_factor = np.where(day_of_week >= 5, 0.5, 1.0)
heat_demand = base_demand * weekend_factor
# Add some random variation
np.random.seed(42)
heat_demand = heat_demand + np.random.normal(0, 5, len(heat_demand))
heat_demand = np.clip(heat_demand, 20, 100)
print(f'Time range: {timesteps[0]} to {timesteps[-1]}')
print(f'Peak demand: {heat_demand.max():.1f} kW')
print(f'Total demand: {heat_demand.sum():.0f} kWh')
# One week, hourly resolution timesteps = pd.date_range('2024-01-15', periods=168, freq='h') # Create realistic office heat demand pattern hours = np.arange(168) hour_of_day = hours % 24 day_of_week = (hours // 24) % 7 # Base demand pattern (kW) base_demand = np.where( (hour_of_day >= 7) & (hour_of_day <= 18), # Office hours 80, # Daytime 30, # Night setback ) # Reduce on weekends (days 5, 6) weekend_factor = np.where(day_of_week >= 5, 0.5, 1.0) heat_demand = base_demand * weekend_factor # Add some random variation np.random.seed(42) heat_demand = heat_demand + np.random.normal(0, 5, len(heat_demand)) heat_demand = np.clip(heat_demand, 20, 100) print(f'Time range: {timesteps[0]} to {timesteps[-1]}') print(f'Peak demand: {heat_demand.max():.1f} kW') print(f'Total demand: {heat_demand.sum():.0f} kWh')
Time range: 2024-01-15 00:00:00 to 2024-01-21 23:00:00 Peak demand: 92.3 kW Total demand: 8007 kWh
In [3]:
Copied!
# Visualize the demand pattern with plotly
demand_series = xr.DataArray(heat_demand, dims=['time'], coords={'time': timesteps}, name='Heat Demand [kW]')
fig = px.line(
x=demand_series.time.values,
y=demand_series.values,
title='Office Heat Demand Profile',
labels={'x': 'Time', 'y': 'kW'},
)
fig
# Visualize the demand pattern with plotly demand_series = xr.DataArray(heat_demand, dims=['time'], coords={'time': timesteps}, name='Heat Demand [kW]') fig = px.line( x=demand_series.time.values, y=demand_series.values, title='Office Heat Demand Profile', labels={'x': 'Time', 'y': 'kW'}, ) fig
Define Gas Prices¶
Gas prices vary with time-of-use tariffs:
In [4]:
Copied!
# Time-of-use gas prices (€/kWh)
gas_price = np.where(
(hour_of_day >= 6) & (hour_of_day <= 22),
0.08, # Peak: 6am-10pm
0.05, # Off-peak: 10pm-6am
)
fig = px.line(x=timesteps, y=gas_price, title='Gas Price [€/kWh]', labels={'x': 'Time', 'y': '€/kWh'})
fig
# Time-of-use gas prices (€/kWh) gas_price = np.where( (hour_of_day >= 6) & (hour_of_day <= 22), 0.08, # Peak: 6am-10pm 0.05, # Off-peak: 10pm-6am ) fig = px.line(x=timesteps, y=gas_price, title='Gas Price [€/kWh]', labels={'x': 'Time', 'y': '€/kWh'}) fig
Build the Energy System¶
The system includes:
- Gas boiler (150 kW thermal capacity)
- Thermal storage tank (500 kWh capacity)
- Office building heat demand
Gas Grid ──► [Gas] ──► Boiler ──► [Heat] ◄──► Storage
│
▼
Office
In [5]:
Copied!
flow_system = fx.FlowSystem(timesteps)
flow_system.add_carriers(
fx.Carrier('gas', '#3498db', 'kW'),
fx.Carrier('heat', '#e74c3c', 'kW'),
)
flow_system.add_elements(
# === Buses ===
fx.Bus('Gas', carrier='gas'),
fx.Bus('Heat', carrier='heat'),
# === Effect ===
fx.Effect('costs', '€', 'Operating Costs', is_standard=True, is_objective=True),
# === Gas Supply with time-varying price ===
fx.Source(
'GasGrid',
outputs=[fx.Flow('Gas', bus='Gas', size=500, effects_per_flow_hour=gas_price)],
),
# === Gas Boiler: 150 kW, 92% efficiency ===
fx.linear_converters.Boiler(
'Boiler',
thermal_efficiency=0.92,
thermal_flow=fx.Flow('Heat', bus='Heat', size=150),
fuel_flow=fx.Flow('Gas', bus='Gas'),
),
# === Thermal Storage: 500 kWh tank ===
fx.Storage(
'ThermalStorage',
capacity_in_flow_hours=500, # 500 kWh capacity
initial_charge_state=250, # Start half-full
minimal_final_charge_state=200, # End with at least 200 kWh
eta_charge=0.98, # 98% charging efficiency
eta_discharge=0.98, # 98% discharging efficiency
relative_loss_per_hour=0.005, # 0.5% heat loss per hour
charging=fx.Flow('Charge', bus='Heat', size=100), # Max 100 kW charging
discharging=fx.Flow('Discharge', bus='Heat', size=100), # Max 100 kW discharging
),
# === Office Heat Demand ===
fx.Sink(
'Office',
inputs=[fx.Flow('Heat', bus='Heat', size=1, fixed_relative_profile=heat_demand)],
),
)
flow_system = fx.FlowSystem(timesteps) flow_system.add_carriers( fx.Carrier('gas', '#3498db', 'kW'), fx.Carrier('heat', '#e74c3c', 'kW'), ) flow_system.add_elements( # === Buses === fx.Bus('Gas', carrier='gas'), fx.Bus('Heat', carrier='heat'), # === Effect === fx.Effect('costs', '€', 'Operating Costs', is_standard=True, is_objective=True), # === Gas Supply with time-varying price === fx.Source( 'GasGrid', outputs=[fx.Flow('Gas', bus='Gas', size=500, effects_per_flow_hour=gas_price)], ), # === Gas Boiler: 150 kW, 92% efficiency === fx.linear_converters.Boiler( 'Boiler', thermal_efficiency=0.92, thermal_flow=fx.Flow('Heat', bus='Heat', size=150), fuel_flow=fx.Flow('Gas', bus='Gas'), ), # === Thermal Storage: 500 kWh tank === fx.Storage( 'ThermalStorage', capacity_in_flow_hours=500, # 500 kWh capacity initial_charge_state=250, # Start half-full minimal_final_charge_state=200, # End with at least 200 kWh eta_charge=0.98, # 98% charging efficiency eta_discharge=0.98, # 98% discharging efficiency relative_loss_per_hour=0.005, # 0.5% heat loss per hour charging=fx.Flow('Charge', bus='Heat', size=100), # Max 100 kW charging discharging=fx.Flow('Discharge', bus='Heat', size=100), # Max 100 kW discharging ), # === Office Heat Demand === fx.Sink( 'Office', inputs=[fx.Flow('Heat', bus='Heat', size=1, fixed_relative_profile=heat_demand)], ), )
Run Optimization¶
In [6]:
Copied!
flow_system.optimize(fx.solvers.HighsSolver(mip_gap=0.01));
flow_system.optimize(fx.solvers.HighsSolver(mip_gap=0.01));
Running HiGHS 1.12.0 (git hash: 755a8e0): Copyright (c) 2025 HiGHS under MIT licence terms
MIP linopy-problem-x9l5ehvr has 2200 rows; 2199 cols; 6740 nonzeros; 336 integer variables (336 binary)
Coefficient ranges:
Matrix [1e-05, 1e+02]
Cost [1e+00, 1e+00]
Bound [1e+00, 5e+02]
RHS [1e+00, 2e+02]
Presolving model
1176 rows, 1008 cols, 2855 nonzeros 0s
840 rows, 672 cols, 3022 nonzeros 0s
840 rows, 672 cols, 3022 nonzeros 0s
Presolve reductions: rows 840(-1360); columns 672(-1527); nonzeros 3022(-3718)
Solving MIP model with:
840 rows
672 cols (336 binary, 0 integer, 0 implied int., 336 continuous, 0 domain fixed)
3022 nonzeros
Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic;
I => Shifting; J => Feasibility jump; L => Sub-MIP; P => Empty MIP; R => Randomized rounding;
S => Solve LP; T => Evaluate node; U => Unbounded; X => User solution; Y => HiGHS solution;
Z => ZI Round; l => Trivial lower; p => Trivial point; u => Trivial upper; z => Trivial zero
Nodes | B&B Tree | Objective Bounds | Dynamic Constraints | Work
Src Proc. InQueue | Leaves Expl. | BestBound BestSol Gap | Cuts InLp Confl. | LpIters Time
0 0 0 0.00% 0 inf inf 0 0 0 0 0.0s
R 0 0 0 0.00% 558.830517 583.6312648 4.25% 0 0 0 464 0.0s
L 0 0 0 0.00% 558.830517 558.830517 0.00% 174 58 0 522 0.1s
1 0 1 100.00% 558.830517 558.830517 0.00% 174 58 0 724 0.1s
Solving report
Model linopy-problem-x9l5ehvr
Status Optimal
Primal bound 558.830516996
Dual bound 558.830516996
Gap 0% (tolerance: 1%)
P-D integral 0.00426289756946
Solution status feasible
558.830516996 (objective)
0 (bound viol.)
1.55431223448e-15 (int. viol.)
0 (row viol.)
Timing 0.13
Max sub-MIP depth 1
Nodes 1
Repair LPs 0
LP iterations 724
0 (strong br.)
58 (separation)
202 (heuristics)
In [7]:
Copied!
flow_system.statistics.plot.balance('Heat')
flow_system.statistics.plot.balance('Heat')
Out[7]:
Storage Charge State¶
Track how the storage level varies over time:
In [8]:
Copied!
flow_system.statistics.plot.balance('ThermalStorage')
flow_system.statistics.plot.balance('ThermalStorage')
Out[8]:
Heatmap Visualization¶
Heatmaps show patterns across hours and days:
In [9]:
Copied!
flow_system.statistics.plot.heatmap('Boiler(Heat)')
flow_system.statistics.plot.heatmap('Boiler(Heat)')
Out[9]:
In [10]:
Copied!
flow_system.statistics.plot.heatmap('ThermalStorage')
flow_system.statistics.plot.heatmap('ThermalStorage')
Out[10]:
Cost Analysis¶
In [11]:
Copied!
total_costs = flow_system.solution['costs'].item()
total_heat = heat_demand.sum()
print(f'Total operating costs: {total_costs:.2f} €')
print(f'Total heat delivered: {total_heat:.0f} kWh')
print(f'Average cost: {total_costs / total_heat * 100:.2f} ct/kWh')
total_costs = flow_system.solution['costs'].item() total_heat = heat_demand.sum() print(f'Total operating costs: {total_costs:.2f} €') print(f'Total heat delivered: {total_heat:.0f} kWh') print(f'Average cost: {total_costs / total_heat * 100:.2f} ct/kWh')
Total operating costs: 558.83 € Total heat delivered: 8007 kWh Average cost: 6.98 ct/kWh
Flow Rates and Charge States¶
Visualize all flow rates and storage charge states:
In [12]:
Copied!
# Plot all flow rates
flow_system.statistics.plot.flows()
# Plot all flow rates flow_system.statistics.plot.flows()
Out[12]:
In [13]:
Copied!
# Plot storage charge states
flow_system.statistics.plot.storage('ThermalStorage')
# Plot storage charge states flow_system.statistics.plot.storage('ThermalStorage')
Out[13]:
Energy Flow Sankey¶
A Sankey diagram visualizes the total energy flows through the system:
In [14]:
Copied!
flow_system.statistics.plot.sankey.flows()
flow_system.statistics.plot.sankey.flows()
Out[14]:
Key Insights¶
The optimization reveals how storage enables load shifting:
- Charge during off-peak: When gas is cheap (night), the boiler runs at higher output to charge the storage
- Discharge during peak: During expensive periods, storage supplements the boiler
- Weekend patterns: Lower demand allows more storage cycling
Summary¶
You learned how to:
- Add Storage components with efficiency and losses
- Use time-varying prices in effects
- Visualize results with heatmaps and balance plots
- Access raw data via statistics.flow_rates and statistics.charge_states
Next Steps¶
- 03-investment-optimization: Optimize storage size
- 04-operational-constraints: Add startup costs and minimum run times